---
title: "Experiment 1"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries}
library(readr)
library(reshape2)
library(plyr)
library(ggplot2)
library(mediation)
library(tidyverse)

```

```{r import}
fulldata <- read_csv("~/Dropbox/Mac/Desktop/Studies/Force/analyses/FINAL/Exp.1/exp1_foranalysis.csv", 
    col_names = FALSE)

```
##Clean data
```{r clean_data}
#Name the columns
names(fulldata)<-c("sub","duration","scenario","pos_neg","force","norm","value")

#Remove rows with any missing data
cleandata<-fulldata[complete.cases(fulldata), ]

#Specify that the variables are factors
cleandata$scenario<-factor(cleandata$scenario,labels=c("ship","doctor"))
cleandata$sub<-factor(cleandata$sub)

#Fix pos_neg values for multiple mediation
clean_numeric<-cleandata
cleandata$pos_neg<-factor(cleandata$pos_neg, labels = c("pos", "neg"))

#separate out the scenarios for multiple mediation
shipmeddata=cleandata[which(cleandata$scenario=="ship"),]

shipmeddata <- shipmeddata %>% mutate(
    norm = as.numeric(norm),
    value = as.numeric(value),
    force = as.numeric(force),
    pos_neg = case_when(
      pos_neg=="pos" ~ 0,
      pos_neg=="neg" ~ 1
    )
)

shipmeddata  <- as.data.frame(shipmeddata)

docmeddata=cleandata[which(cleandata$scenario=="doctor"),]
docmeddata <- docmeddata[-5,] ## dunno why but row five has (subj 10) has an NaN for norm..

docmeddata <- docmeddata %>% mutate(
    norm = as.numeric(norm),
    value = as.numeric(value),
    force = as.numeric(force),
    pos_neg = case_when(
      pos_neg=="pos" ~ 0,
      pos_neg=="neg" ~ 1
    )
)

docmeddata  <- as.data.frame(docmeddata)

fullmeddata=cleandata
fullmeddata<-fullmeddata %>% mutate(
  norm = as.numeric(norm),
  value = as.numeric(value),
  force = as.numeric(force),
  pos_neg = case_when(
    pos_neg=="pos" ~ 0,
    pos_neg=="neg" ~ 1
  )
)
fullmeddata <- as.data.frame(fullmeddata)


#separate out the scenarios
shipdata=cleandata[which(cleandata$scenario=="ship"),]
docdata=cleandata[which(cleandata$scenario=="doctor"),]
```

##Descriptives
```{r descriptives}

#get descriptives
across_scen_desc <- ddply(cleandata, c("pos_neg"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(across_scen_desc)

force_desc <- ddply(cleandata, c("scenario","pos_neg"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc)

norm_desc <- ddply(cleandata, c("scenario","pos_neg"), summarise,
               N    = length(norm),
               mean = mean(norm),
               sd   = sd(norm),
               se   = sd / sqrt(N)
)
print(norm_desc)

value_desc <- ddply(cleandata, c("scenario","pos_neg"), summarise,
               N    = length(value),
               mean = mean(value),
               sd   = sd(value),
               se   = sd / sqrt(N)
)
print(value_desc)

```

##Final plots
```{r plots}

###FINAL BAR CHART####
#Create a new ordered condition vairable for chart labels
cleandata$Condition<-factor(cleandata$pos_neg, ordered=TRUE, levels=c("neg","pos"), labels = c("Morally Wrong", "Morally Neutral"))

exp1_fig1a<-cleandata %>%
  ggplot(aes(Condition,force, fill=Condition)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2, position='dodge') +
  geom_point(aes(Condition, color=Condition), size=2,
             position=position_jitterdodge(dodge.width=.9, jitter.width = .8), alpha = .5) +
  scale_color_manual(values = c("Morally Neutral" = "#017073", "Morally Wrong" = "#8a433e")) +
  scale_fill_manual(values=c("Morally Neutral" = "#00BFC4", "Morally Wrong" = "#F8766D")) +
  xlab("Moral Valence of the Actual Action") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 75)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=20, vjust=-1),
        axis.title.y = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "right")

#ggsave("exp1_fig1a.jpg",height = 10, width= 7, dpi=600)

###FINAL SCATTER PLOTS###
exp1_fig1b<-ggplot(cleandata,aes(y=force,x=value))+
  geom_point()+
  geom_smooth(method=lm) + 
  xlab(str_wrap("How strongly do you agree that it would have been A GOOD IDEA to do something else instead?", 50)) + 
  scale_x_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 30)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=16, vjust=-1),
        axis.title.y = element_text(size=16, vjust=.35),
        panel.background = element_blank(),
        axis.line = element_line(colour = "grey"),
        #panel.grid = element_blank(),
        legend.position = "right")

#ggsave("exp1_fig1b.jpg",height = 5, width= 8, dpi=600)

exp1_fig1c<-ggplot(cleandata,aes(y=force,x=norm))+
  geom_point()+
  geom_smooth(method=lm) + 
  xlab(str_wrap("How strongly do you agree that it would have been UNUSUAL to do something else instead?", 50)) + 
  scale_x_continuous(breaks=c(0, 25, 50, 75, 100), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 30)) +
  scale_y_continuous(breaks=c(0, 25, 50, 75, 100), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title.x = element_text(size=16, vjust=-1),
        axis.title.y = element_text(size=16, vjust=.35),
        panel.background = element_blank(),
        axis.line = element_line(colour = "grey"),
        #panel.grid = element_blank(),
        legend.position = "none")

#ggsave("exp1_fig1c.jpg",height = 5, width= 8, dpi=600)

```
##Final Regressions
```{r testing effect of condition echo=FALSE}
#Subject level analysis: Predict force judgments by condition, with a random intercept for scenario. 
model<-lme4::lmer(force~pos_neg + (1|scenario), data =cleandata)
summary(model)

#Drop condition from the model and compare to the original model
smmod<-lme4::lmer(force~(1|scenario), data =cleandata)
anova(model,smmod)
```

```{r testing effect of unusualness and value echo=FALSE}
#Subject level analysis: Predict force judgments with norm and value judgments with a random intercept for scenario

library(car)
mod1a<-lme4::lmer(force ~ norm + value + (1|scenario), data=cleandata)
summary(mod1a)

#Test the effect of unusualness by dropping norm from model and comparing it to original model
mod1b<-lme4::lmer(force ~ value + (1|scenario), data=cleandata)
summary(mod1b)
anova(mod1a, mod1b)

#Test the effect of value from dropping value from model and comparing it to the original model
mod1c<-lme4::lmer(force ~ norm + (1|scenario), data=cleandata)
summary(mod1c)
anova(mod1a, mod1c)

#To test the difference in strength of the absolute value of the betas, create a new variable for value that is reverse coded. 
cleandata<-mutate(cleandata,value_rev=100-value)
mod2<-lm(force ~ norm + value_rev + scenario, data=cleandata)
summary(mod2)

#run Wald test to compare relative strength of betas from original model
linearHypothesis(mod2, c("norm = value_rev"))
```

##Final mediation analyses

```{r full multiple mediation echo=TRUE}

#Run multiple mediation testing the mediating effect of unusualness, controlling for value, on the relationship between condition and force judgments (with scenario as a covariate)
multi_med_normmain<-multimed(outcome="force", med.main="norm", med.alt="value", treat="pos_neg", covariates = "scenario", data=fullmeddata, sims = 1000)
summary(multi_med_normmain)

#Run multiple mediation testing the mediating effect of value, controlling for unusualness, on the relationship between condition and force judgments (with scenario as a covariate)
multi_med_valmain<-multimed(outcome="force", med.main="value", med.alt="norm", treat="pos_neg", covariates = "scenario", data=fullmeddata, sims = 1000)
summary(multi_med_valmain)

```
